home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
rwvector.lha
/
RWVector2.1
/
src
/
dmath.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1989-08-18
|
3KB
|
141 lines
/*
* Double Precision math function definitions
*
* Copyright (C) 1988, 1989.
*
* Dr. Thomas Keffer
* Rogue Wave Associates
* P.O. Box 85341
* Seattle WA 98145-1341
*
* Permission to use, copy, modify, and distribute this
* software and its documentation for any purpose and
* without fee is hereby granted, provided that the
* above copyright notice appear in all copies and that
* both that copyright notice and this permission notice
* appear in supporting documentation.
*
* This software is provided "as is" without any
* expressed or implied warranty.
*
*
* @(#)dmath.cc 2.1 8/18/89
*/
#include "rw/DComplexVec.h"
#include <math.h>
DComplexVec
rootsOfOne(int N, unsigned nterms)
{
/* Use the recursion formulae
* sin(a+b) = sin(a)cos(b) + cos(a)sin(b)
* cos(a+b) = cos(a)cos(b) - sin(a)sin(b)
*/
if(!N){
RWnote("::rootsOfOne()", "Zero'th root requested.");
RWerror(DEFAULT);
}
unsigned order = abs(N);
double b = 2.0*M_PI/N;
register double cosb = cos(b);
register double sinb = sin(b);
register double cosa = 1.0; // implies a=0
register double sina = 0.0;
DComplexVec roots(nterms);
if(nterms){
roots(0) = DComplex(cosa, sina);
register double temp;
for(register int i = 1; i<nterms; i++){
temp = cosa*cosb - sina*sinb;
sina = sina*cosb + cosa*sinb;
cosa = temp;
roots(i) = DComplex(cosa, sina);
}
}
return roots;
}
/*
* Expand a complex conjugate even series into its full series.
*/
DComplexVec
expandConjugateEven(const DComplexVec& v)
{
unsigned N = v.length()-1;
DComplexVec full(2*N);
// Lower half:
full.slice(0,N+1,1) = v;
// Upper half:
full.slice(N+1,N-1,1) = conj(v.slice(N-1,N-1,-1));
return full;
}
/*
* Expand a complex conjugate odd series into its full series.
*/
DComplexVec
expandConjugateOdd(const DComplexVec& v)
{
unsigned N = v.length()-1;
DComplexVec full(2*N);
// Lower half:
full.slice(0,N+1,1) = v;
// Upper half:
full.slice(N+1,N-1,1) = conj(-v.slice(N-1,N-1,-1));
return full;
}
/*
* Expand an even series into its full series.
*/
DoubleVec
expandEven(const DoubleVec& v)
{
unsigned N = v.length()-1;
DoubleVec full(2*N);
// Lower half:
full.slice(0,N+1,1) = v;
// Upper half:
full.slice(N+1,N-1,1) = v.slice(N-1,N-1,-1);
return full;
}
/*
* Expand an odd series into its full series.
*/
DoubleVec
expandOdd(const DoubleVec& v)
{
unsigned N = v.length()+1;
DoubleVec full(2*N);
// Lower half:
full.slice(1,N-1,1) = v;
// Upper half:
full.slice(N+1,N-1,1) = v.slice(N-2,N-1,-1);
full(0) = 0;
full(N) = 0;
return full;
}